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down on computer time, a search was made up over only the two variables 
v,w, since u, the heading angle, is of relatively small significance in 
the computation of the Hamiltonian. The search consisted of taking all 
possible combinations of v,w where v ranged over the set of values 3,6, 
9,12,15 and w ranged over the values 200,400,600,800,1000 and comparing 
the Hamiltonian for each set with the one generated in the numerical rou- 
tine. It is possible for the search with this size grid to fail to de- 
tect a set of controls that yield a smaller value for the Hamiltonian. 
However, the time factor is critical and this grid size was felt to fur- 
nish a satisfactory compromise. A search relying on gross computation 
can easily become completely unrealistic in terms of the computer time 
required. 

The Hamiltonian 



(3.9) H = A • V 

has a constant "^alue on an extremal, with allowance being made for round 
off errors, whenever t does not occur explicitly in H. 

At this point, we might mention that there are two types of 
tions of the control variables in the classical literature, called weak 
variations and strong variations. [4 *] Weak variations are variations in 
which the | c/u^ | are ’’small’’ for each time stej^^strong variations are 
variations in which |</u^| dt is ”small^\ That is, in weak variations 
only values of control near those used are compared but if strong varia- 
tions are considered, then the new control function may not be ’’near” the 
one used. 

All methods of determining the routes are methods of variations which 
deform a given path. A path which furnishes a relative minimum, if we 
allow only we variations , may not furnish such a minimum if strong var- 
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iations are allowed, as will be seen in section 6, 

5. The Numerical Routine. 

The routine for determining the route is given in this section. 
Heuristic Discussion : 

Let us guess a set of values for the parameters h^, h^* We will 
then use this set of values to determine the control variables for each 
time t by the minimum principle to determine a route. The terminal point 
thus generated will, in general, differ from the desired one. By chang- 
ing the values of h^, h^ appropriately, this terminal point will be 
forced toward the desired point x^, y^. 

Mathematical Derivation ; 

First we see from (4.3) that 

(5.1) SA ^ A^dh^ + 

Note that "'/ = 1, and this in turn implies <TV = 0. These facts will be 
used in the following equations. 

Next, from the first of the Euler equations (4.4) by taking the 
total differential we find that 

- V sin u/A + V cos u//” + (-Av cos u -.^v sin u + f ) 

(5.2) 

+(-A sin u - y cos u f )^v + f <Tw = 0, 

uv uw 

if we assume that we can change A, , u, v, w, at fixed x, y, z, t. But 
since 

H = - A V cos u V sin u + f 
uu uu 

H =-A sinu-/^ cos u + f 
uv uv 

H = f 
uw uw 

(5.2) becomes 
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ABSTRACT 



In this paper the details of computing an optimum route for a sub- 
marine are studied. 

Typical functions representing the listening devices were used. It 
was found that in some cases several extremals existed and it was neces- 
sary to set up tests for the Legendre and Vfeierstrass conditions. The 
problem is further complicated by the fact that the optimum control- 
variables may lie on the boundary of the region of allowed values and 
further routines must be adjoined for this. Further, corners may occur 
and in particular the control may move discontinuously from a boundary 
point to an interior point or vice versa. The routines were made up to 
effect a compromise between the need for accuracy and reasonable com- 
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Introduction 



The purpose of this paper was to develop a numerical routine for 
solving the submarine routing problem; the problem is for the submarine 
to choose a course that minimizes the probability that it will be de- 
tected. Several difficulties arise in the numerical solution and sub- 
routines were included to take care of these. 

In type, the problem is a problem of Bolza in calculus of varia- 
tions. It is complicated by the fact that there may be several routes 
each of which appear to be the solution. They all begin and end at the 
desired points and they all satisfy the Euler equations. It is only 
after routines are incorporated to check other conditions, the condi- 
tions of Legendre and Weierstrass, that it becomes clear whether a 
particular solution is the desired solution. When these conditions are 
not satisfied this fact must be determined and a routine made up to 
determine the controls to satisfy it. 

In general, the route is generated as an extremal, a solution to 
the Euler equations, though the basic principle is that the control must 
minimize the Hamiltonian. The problem is complicated by the fact that 
the control which furnishes the minimxjm may lie on the boundary to the 
region, as when the submarine is at maximum depth. Several subroutines 
must be adjoined to treat the case when the control lies on one of the 
bounding faces or edges. 

The existence of a corner introduces further complications; at a 
corner the control is a discontinuous function of time. It is necessary 
to make up a search routine for other values of the control which may 
decrease the Hamiltonian, and it is necessary to compromise between the 
demands of computing time and accuracy in this routine. 
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!• Statement of the Problem. 



The basic problem studied here is the following. A submarine 
located at point is to make a voyage to point x^,y^ in a specified 

time T. Throughout this voyage the submarine is subjected to enemy de- 
tection devices whose capabilities are assumed known statistically. If 
the submarine has previously gone undetected, let the probability of 
detection in a time interval At be approximately f (K,y ,u,v,w, t) At, and 
hence the probability of being detected along the route satisfies the 
equation 

(1.1) dp = (1 - p) f (x,y,u,v,w,t) dt. 

The function f is the best estimate of the enemy's detection capabili- 
ties based on information we have about his listening devices, tests we 
have run on our submarines using comparable devices, the distance in- 
volved in the trip, and other information available to us. 

Equation (1.1) can be simplified by letting 

(1.2) z = - ln(l - p) 
which gives 

(1.3) z = f. 

Since we are primarily interested in long routes, and the time to 
change depth, speed, and heading angle is assumed small compared to 
total time, it will be ignored. 

2. Equations. 

Because routes of approximately 2500 miles are of primary interest, 
the flat earth assumption will be used. The equations governing the 
system may then be written as 

X = V cos u 
(2.1) 7 = V sin u 

z = f (x,y ,u,v,w,t). 
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limited range. The constant reflects the possibility of detection by 
some passing ship, for example. 

The thermocline is defined as the depth at which the temperature 
gradient (rate of decrease of temperature with increasing depth) is a 
maximxjm. In equations (2.3), Wq is the depth of the thermocline, s^ is 
the distance from the terminal point (x^,y^) of the route to the main 
concentration of enemy listening devices which is represented by the 
point (X 2 ,y 2 )> and 0 is the angle that the perpendicular to the enemy 
shore line makes with the x-axis. The constants a^,b^,Qj,dj,and Wq are 
chosen so as to simulate, by the function f, conditions as they exist 
in any given situation. 

The two functions f^ and f^ are intended to be typical of the 
functions which do describe the enemy’s defenses. The function f^, 
representing the passive enemy defense, tends to decrease as the dis- 
tance from the enemy shoreline increases. This mathematical model of 
the enemy’s passive defenses also tends to be insensitive in the region 
of the thermocline. The function which represents the enemy active 
defense, was constructed so as to emphasize the enemy’s surface search. 
For this reason, this function decreases as the depth increases. 

If submarine routing as described in this paper were to be made a 
part of naval operations, generating the function f would be a problem 
for ihteiligence and engineers. Such things as the sea state and the 
nature of the ocean floor in the region where the submarine is operating 
would then have to be included in the function f. These functions would 
also vary with time, reflecting changing sea state, etc. 

The problem now becomes that of determining the control variables 
as functions of time to effect the desired optimization. That is, we 
want to choose the heading, the speed, and the depth so as to go from 
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the initial point (x^jy^) to the terminal point (x^,y^) with the minimum 
probability of detection. 

By (1.1) the probability of being detected along the route is 

(2.5) p(T) = 1 - exp 1^- z(T) J 
and hence this is the quantity we wish to minimize. 

3. Adjoint Equations. 

Let us consider any route and a neighboring route. The neighbor- 
ing route will be generated by replacing u,v,w on the original route 
by u + J*u, V +cTv, w + <Pw respectively. This generates first-order 
changes in x,y,z satisfying the differential equations 



cTx = - V sin u cTu + cos u cTv 

(3.1) cTy = V cos u cTu + sin u <fv 

cTz = f cTu + f </'v + f o^w + f cTx + f <Ty. 
u V w X y 



a f 

* 3u 



The notation f„ ^ 9 etc. as used above will be employed throughout 



the remainder of the paper. 

Let us now introduce three Lagrange multipliers which are 

some unspecified function of time t. Multiply equations (3.1) by A , , 
V in that order, add the three equations, and integrate the result over 
the interval (0,T). These operations yield 
T r 

(1 r + V sin u cTu - cos u a'v +/'(<Ty - 

Jf L 



V cos u ^u - sin u «^v) + V(rf'z - f /x - f (fy 

X y 

-f (Tu-f<A^-f ^w) dt. 

u V w J 

Separating the terms containing the variations of the state variables 
from those containing the variations of the control variables, we get 
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(3.3) 



^ <Tx + <Ty + v'(<Tz 



f «Tx - 

X 




dt 



[- 



^vsinu+Z^vcos u+Vf 



u]'^" * [ 



}[ cos u 



+ /'sin u + Vf 1<Tv + Vf </*w) dt. 
vj w 

Integrating by parts in (3.3) just the terms on the left involving 
derivatives of state variables, we get 



(3.4) 



<(x +J'Jy +y<fz~^ - ^ ^erx(^ 



+ Vf ) 
X 



+ (.^ + vf ) + <Tz''^J 



dt 



which combined with the right-hand equations of (3.3) yields 

T ^ 

^ ^<^x Jy <(zj - \ <Ac(3v + f^) + 



(3.5) (/y(J^ 



0 '0 



^ V sin u +>^v cos u 



+Vf ) /u + ( ^ cos u + >^sin u + Vf ) <Tv + Vf (Tw dt . 
u V w J 

Now, choose A V in such a manner that they satisfy the differential 

equat ions 

^ = -Vf 

X 

(3.6) ^ = -*Vf 

- y 

= 0 . 



Equations (3.6) are called the adjoint equations of the variation- 
al equations (3.1). With this choice of A ,V equations (3.1) reduce 
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to 



(3.7) 



T 

/i (Tx + >^ «Ty + V <Tzj = \ V sin u - /v cos u 

0 -^0 

r 

-Vf^ j /u dt + \ cos u+)^sinu + V /v dt 

-'o 



-^0 



<Tw dt. 



Note that this formula gives us a relation between the terminal values 
of x,y,z and an integral made up of the variations of the control varia- 
bles cTu, cTv, (i\j. Note also the important fact that the values of 
cTx, </y, iPz interior to the interval (0,T) are not needed. 

For convenience in the sections to follow, the vector notation 




will be used. The scalar product of these two vectors 



(3.9) H = A * V 
is called the Hamiltonian. 

Also for future use, let us define the following three particular 
solutions to the adjoint equations 
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4. Conditions for a Minimum. 

In this section, the necessary condition for a minimum will be 
given. 

If a path is to provide the desired minimum, it must first be admis- 
sible. 

Admissibility . One requirement for a route to be admissible is that 
it begin and end at the desired points, i.e., x(0) = 0, y(0) = 0 and 
x(T) = x^, y(T) = y^ for some solution to the differential equations 



X = V cos u 
(2.1) y = V sin u 
z = f . 

A further condition for admissibility is that the depth and speed satisfy 
the inequalities 

0 f; w i w 



(4.1) 



max 



0 6 . V ^ V 



max 



where w and v depend upon the specific class of submarine under 
max max 

consideration. 



A set of control variables which are piecewise continuous and sat- 
isfy (4.1) are called allowable. Note that allowability is a local con- 
straint, in terms of time, on the set of control variables. A path which 
satisfies the above constraints, (2.1), (4.1), is admissible . 
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Our problem is to find the one route, among all admissible routes, 
such that z(T) is a minimum. Minimizing z(T) in turn minimizes p(T) and 
this is our objective. 

For a path to furnish a minimum, it must be admissible and also sat- 
isfy the following conditions: 

1. Euler Equations. 

2. Legendre Condition. 

3. Weierstrass Condition. 

4. Envelope Condition. 

The envelope condition was not investigated and will not be treated in 
this paper. The order of the conditions as given above is used since 
this is the usual order in which they will be checked in a problem. 

A point that should be kept in mind throughout is that on a path 
which furnishes the desired minimum, the Hamiltonian must be minimized 
at each value of t. This is the Weierstrass-Pontryagin maximum princi- 
ple with a change of sign. To satisfy this one criterion, the Euler 
equations, the Legendre condition, and the Weierstrass condition must all 
be satisfied. 

The Euler equations, a condition on the first derivatives, require 
that the Hamiltonian have a stationary value. 

In the Legendre condition the second derivatives are checked for a 
minimum. [ 2 J The Euler equations and the Legendre conditions are both 
local conditions on the control variables u,v,w. 

Finally, in the Weierstrass condition we compare the Hamiltonian 
for the values of u,v,w used with the Hamiltonian for all allowable values 
of u,v,w, for all values of t. 

Euler Equations . For a minimum, the control variables must be chosen 
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so they minimize 

(3.9) H = A * V 

as compared with all allowable controls, for all t, 0<t^T, for some 
solution to the adjoint system of equations. Let us consider our solu- 
tion A to the adjoint equations in the form 



( 4 . 2 ) A = A j + 1^2 X2 + 4 ^ A3 

where h^, h2> h^ are arbitrary constants and A^, A3 are defined 

in equations ( 3 . 10 ). 

We are now faced with the problem of having three constants in our 
solution A to the adjoint equations. But we may choose one relation 
among the constants h^, h^, h^. We chose h^ to be 1, and then 

( 4 . 3 ) A = A J + h2 A2 + A3. 



This then leaves us with the problem of finding the other two constants 
so that x(T) , y(T) will assume the desired terminal values x^, y^. 

If the values of the control variables lie inside the domain of 
allowed values, then the Euler equations 



-^-S— = H = -^ V sin u +/'v cos u +i^f = 0 

ju u u 



dH 



(4.4) " " = H = ^cos u + /' sin u + Vf = 0 

3 V V V 

-^-2- = H = Vf =0 
3 w w w 

are the first necessary conditions for the desired minimum. The Euler 
equations when combined with the adjoint equations (3.6) are referred 
to as the Euler-Lagrange equations. Solutions to the equations of motion 
which also satisfy the Euler-Lagrange equations are called extremals . 

The Euler equations, however, do not insure that we effect the desired 
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minimum for H; they are only necessary conditions. It is then necessary 
to investigate additional conditions, the next being the Legendre condi- 
tion. 



Legendre Condition . This is an investigation of the second-degree 

terras of the Taylor expansion of the Haaiiltonian. If we can expand H in 

a series in du, dv, dw, valid in some neighborhood of (0,0,0), the first 

necessary condition for H to be a minimum is that H = H = H = 0, as 

u V w 

stated in equations (4.4). Ihe second necessary condition is that the 
quadratic form 







/H 


‘h 


H \ 


^ du \ 






/ uu 




uw \ 


(4.5) 


(du dv dw) 


H 


H 


H 


dv 


vu 


w 


vw 








\H 


H 


H / 


\ 1 






\ wu 


wv 


ww ] 


\ / 



be positive semi-definite at least; hope it will be positive definite. 
This condition, which is known as the Legendre condition, [3^ may be ex- 
pressed 



(4.6) 



H 


> 0 




uu 






H 


H 




uu 


uv 




H 


H 


> 


vu 


w 






H, 


H 


uu 


uv 


uw 


H 


H 


H 


vu 


w 


vw 


H 


H 


H 


wu 


wv 


ww 



> 0 



for the case in which the control variables u,v,w all lie interior to 
the domain of allowed values. 

If one or more of the control variables are on the boundary, the 
conditions are altered accordingly. Consider the case in which w = w 



max 



for some part of the route. The conditions then become 
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H < 0 



(4.7) 



w 




H 


H 


uu 


uv 


H 


H 


vu 


w 



> 0 



whenever w = If w = and = 0, then we must check to ensure 

that >0. Of course if H^> 0 the minimum is not on the boundary at 
that point. 

If V = V , conditions equivalent to those above will be used, 
max 

When both w = w and v = v , the conditions read 
max max 



H < 0 

w 

(4.8) H 0 

V 

H > 0, 
uu ’ 

and if H = 0 on the boundary then we must check H >0, and similarly 
w ww 

for H and H . 

V w 

If the Euler equations are satisfied and if equations (4.6) are sat- 
isfied at some interior point of u,v,w space, then these values yield a 
local minimum; H is smaller at that point than at any other point u,v,w 
in Some neighborhood. However, the point may not furnish the minimum 
value for H; it is necessary in theory to compare the value with the value 
for all allowable values of u,v,w. This condition, that the value chosen 
minimize H, is called the Weierstrass condition . 

A comparison must be made between the Hamiltonian for the u,v,w gen- 
erated, H(u,v,w), and for all other allowable combinations of control 
variables, say, u^, v^, w^. If H(u,v,w) > H(u^, v^, w^) , then the 
controT variables u,v,w must be replaced by the new set u^, v^, w^, 
which yield the smaller value for the Hamiltonian. 

Unfortunately there is no easy way to check this condition. To cut 
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down on computer time, a search was made up over only the two variables 
v,w, since u, the heading angle, is of relatively small significance in 
the computation of the Hamiltonian. The search consisted of taking all 
possible combinations of v,w where v ranged over the set of values 3,6, 
9,12,15 and w ranged over the values 200,400,600,800,1000 and comparing 
the Hamiltonian for each set with the one generated in the numerical rou- 
tine. It is possible for the search with this size grid to fail to de- 
tect a set of controls that yield a smaller value for the Hamiltonian. 
However, the time factor is critical and this grid size was felt to fur- 
nish a satisfactory compromise. A search relying on gross computation 
can easily become completely unrealistic in terms of the computer time 
required. 

The Hamiltonian 



(3.9) H = A • V 

has a constant ^lalue on an extremal, with allowance being made for round 
off errors, whenever t does not occur explicitly in H. 

At this point, we might mention that there are two types of 
tions of the control variables in the classical literature, called weak 
variations and strong variations. ^4 ^ Weak variations are variations in 
which the ] | are ’’small” for each time stejj^^strong variations are 

_T * 

variations in which \ dt is ’’small”. That is, in weak variations 

-'O 

only values of control near those used are compared but if strong varia- 
tions are considered, then the new control function may not be ’’near” the 
one used. 

All methods of determining the routes are methods of variations which 
deform a given path. A path which furnishes a relative minimum, if we 
allow only wel^ variations , may not furnish such a minimum if strong var- 
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iations are allowed, as will be seen in section 6. 

5. The Numerical Routine. 

The routine for determining the route is given in this section. 
Heuristic Discussion ; 

Let us guess a set of values for the parameters h^, h 2 * We will 
then use this set of values to determine the control variables for each 
time t by the minimum principle to determine a route. The terminal point 
thus generated will, in general, differ from the desired one. By chang- 
ing the values of h^, h^ appropriately, this terminal point will be 
forced toward the desired point y^. 

Mathematical Derivation ; 

First we see from (4.3) that 

(5.1) SA = A^dh^ + 

Note that i, and this in turn implies ~ 0. These facts will be 

used in the following equations. 

Next, from the first of the Euler equations (4.4) by taking the 
total differential we find that 

- V sin u/^ + V cos ucT/' + (-Av cos u -/'v sin u + f ) (Tu 

(5.2) 

+(-A sin u - cos u f )<Tv + f <Tw = 0, 

uv uw 

if we assume that we can change A, , u, v, w, at fixed x, y, z, t. But 

a: 

since 

H = V cos u -Z' V sin u + f 
uu uu 

H = - sin u - ^ cos u + f 
uv uv 

H = f 
uw uw 

(5.2) becomes 
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(5.3) 



(5.3) 


-V sin u + V cos u J'/' + H <Tu + H + H A; = 0. 

uu uv uw 


Similarly 

(5.4) 


cos u + sin u J'/' + sin u +/'cos u + f ) /u 

vu 


and since 


+ f rf'v + f = 0 

w vw 




H =-P^sinu+^cosu + f 

vu vu 




H = f 

W VV 


(5.4) becomes 


H = f 

vw vw 


(5.5) 


cos u/A + sin u <f}' + H c/u + H cTv + H <Tw = 0. 

vu VV vw 



The third Euler equation gives us 



(5.6) 
But since 


f cTu + f <^v + f cTw = 0. 
wu wv ww 




f = H 

wu wu 




f = H 

wv wv 




f = H 

ww ww 



(5.6) becomes 



(5.7) 


H <Tu + H <A^ + H <fw = 0. 

WU WV WW 



Now consider the equations (5.3), (5.5), (5.7) as three equations in the 
three unknowns cTu, <fv, 





H </^u + H <A/+H ^w = vsinu/^-vcosu</V 
uu uv uw 


(5.8) 


H <Ai + H <A/ + H <fw = -cos u - sin u 

vu w vw 




H cAi + H ^v+H <Tw=0. 
wu wv ww 
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If the determinant of coefficients of cTu, dvy Xw in (5.8) does not vanish, 
we can solve for </u, <Tv, </w by using Cramer’s Rule as follows: 

Denote the determinant of coefficients by D, i.e., 



(5.9) D= 



Then 



Su = 



(5.10) 



«uu «uv «uw 

H H H 
vu w vw 

H H H 
WU WV WW 



(v sin u - V cos u </V) 

(- cos u - sin u «/V) 

0 



H H 
uv uw 

H H 
w vw 

H H 

WV WW 



D 




Notice that <Tw could be found in the same manner, but is not needed 

since it is not used in the numerical routine. Next, since from (5.1) 

= dh, and //' = dh 
1 2 






(5.11) 



(Tu = -(v sin u dh, - v cos u dh.) (H H - H H ) 
1 2 w WW WV vw 

D 

- (cos u dh- + sin u dh„) (H H - H H , ) 

1 2 uv WW WV uv 

D 
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/v = (v sin u dh^ - V cos u dh^) 

D 

+ (cos u dh + sin u dh.) (H' „ H - H H ) 
1 2 uu ww wu uw 

D 

To simplify the form of the equation let us set 



= f(H H - H H ) (- V sin u) + (H H 

L w v/w wv vw uv ww 

H H ) (- cos u)l / D 
wv uw J 

12 r 

S = (H H - H H ) (v cos u) + (H H 

L w ww wv vw uv ww 

(5.12) H H ) (- sin u)l / D 

wv uw J 

21 " 

S = 1 (H H - H H ) (v sin u) + (H,„ H 

L vu ww wu vw uu ww 

H H ) (cos u)1 / D 
wu uw J 



= f(H H - H H ) (- V cos u) + (H H 
L vu ww wu vw uu ^ 

H H ) (sin u)l / D; 
wu uw J 



then the above equations may be rewritten 



(5.13) 



r oil ol2 

Ju = S dh^ + S dh^ 

/. 21 22 

dv = S dh- + S dh„ 
1 2 



We now have the machinery to derive </’x(T) , JyiT). In deriving 
(fxCT) , we use equations (3.5) and a particular choice for A » namely 
the Aj ~ \ I defined in equations (3.10). This choice for A in 



equation (3.7) yields 




<Tv dt - 




V sin u 



(Tu dt 
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which becomes, after noting that (T^CO) = 0 and substituting <fu, <^v 
from equations (5.13), 

r-T 



r 

= \ j^cos u 



(5.14) cTx(T) = \ [cos u (S^^ dhj^ + dh^) - v sin u (S^^ dh^ 



+ dh^)] dt. 



0 



Similarly, using A o ~1 ^ / > find 



(5.15) <Ty(T) 



r 

= I ^sin u 

-'o 



21 22 11 
(S dh + S dh ) + V cos u (S dh 



12 1 
+ S dh2)Jdt. 



Equations (5.14) and (5.15) can be rewritten as 



■■r 



21 11 
(Tx(T) = dhj^ I ^cos uS -vsinuS Jdt + 

0 



dh. 



(5.16) 



r r • ol2 1 

\ cos u S - V sin u S I 



(fy(T) = dh 



1 [“" " 

n 



21 

S + V cos 



11 T 
u S J 



dt 



dt + 



r 



dh^ j jsin u 
0 ^ 
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S + 



12 -] 
V cos U S I 



dt. 



If we set 



■i; 



21 11 

= \ (cos u S - V sin u S ) dt 



■t 



22 12 , 

Ai 2 = \ (cos uS -vsinuS )dt 
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21 



"22 



f 



21 11 
(sin u S + V cos u S ) dt 






22 12 
(sin u S + V cos u S ) dt , 



equations (5.16) become 

(Tx(T) = dh^ + ^^^2 

(5.17) 

cfy(T) = dh^ + dh^. 

Equations (5.17) give us the mechanism for a Newton - Raphson iteration 
scheme for correcting h^, h 2 , if we make the following substitutions 

cTx(T) = - x(T) 

(5.18) 

cfy(T) = - y(T) . 

Note that in equations (5.18) (Tx(T) , <Ty(T) were equated to the desired 
terminal values minus those which were attained. Using equations (5.18), 
equations (5.17) become 

Xrj, - x(T) - \i + A ^2 ^^2 

(5.19) 

yrj, - y(T) = dh^ + A^^ ^^2 

from which we are able to generate corrections for the constants h^, h 2 . 

Note above that in equations (5.19) the coefficients A^ 2 > ^21’ 

A 22 are integrals with respect to time. In the numerical routine, this 
integration is accomplished by a Runge-Kutta numerical integration scheme, 
contained within which is a Newton-Raphson iteration to generate the neces- 
sary changes in u,v,w to accomplish the integration. 

The mathematical basis for the Newton-Raphson iteration scheme to 
generate corrections for u, v, w is the following. From the equations 



H 


du 


+ H 


dv 


+ H 


dw = - H 


uu 




uv 




uw 


u 


H 


du 


+ H 


dv 


+ H 


dw = - H 


vu 




vv 




vw 


V 
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H , du + H dv + H dw 
wu wv ww 



= - H 



w 



we find 



du = 



- H 


H 


H 


u 


uv 


uw 


- H 


H 


H 


V 


w 


vw 


- «w 


Hwv 





(5.21) 



where 









- «u 


Huw 






\u 


- H 

V 


«vw 


dv = 






- «w 


Hww 








D 








H 

uu 


H - H 

uv u 






H 

vu 


H - H 

W V 


dw = 




H 

wu 


H 

wv 


• H 

w 








D 








«uu 


\v 


H 

uw 


D = 




«vu 


H 

w 








Hwu 


^wv 




scheme 


then has the 


form 




= ^n-l 


+ du^ 






3 ^ 


'^n-l 


+ dv^ 






1 — i 

1 

II 


+ dw 

n 





(5.22) 



where n is the index for the Newton-Raphson iteration; du, dv, dw are 
those found in equations (5.21). 
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6. Control Variables on the Boundary. 

In section 4, it was noted that for a path to be admissible the con- 
trol variables v, w had to satisfy the inequalities 





o 

l^ 


w w 


(4.1) 




max 




0 ^ 


V — V 

max 



Let us consider the situation in which the depth w assumes the max- 
imum depth w for part or all of the route and the other bounded con- 
trol variable v remains within its prescribed bounds. We must amend the 
routine for determining the controls as follows. We take the maximum 
depth of the type of submarine being considered and read this infor- 

mation into the program of the numerical routine which calculates the 
route. If the Newton-Raphson iterations as established in (5.22) yield 
a value for w which exceeds w , we set w equal to w and generate 
u and V, using the following iteration scheme. From our Newton-Raphson 
iteration equations, 



( 6 . 1 ) 



H du + H dv = - H 
uu uv u 



H du + H dv = - H , 

VU W V 



we get corrections to the controls, 



( 6 . 2 ) 



du = 



dv 



- H H 
u uv 

-«v «w 



H 

uu 


H 

uv 


H 

VU 


H 

w 


H 

uu 


- H 

u 


H 

VU 


- H 

V 




H 

uv 


\u 


H 

vv 
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Then the iterations discussed in the previous section take the form 



(6.3) 



u = u , + du 
n n-1 n 



= V T + dv„ . 

n n-1 n 



Thus we may modify our Newton-Raphson iteration scheme for the con- 
trols when the minimum value of H occurs for w on the boundary. The 
modified Newton-Raphson iterations generate successive corrections to u 
and V so as to produce an admissible path. The numerical routine is pro- 
grammed so that it is possible for the depth to assume its maximum for 
some part of the route without remaining fixed at maximum depth after 
once assuming it. In section 7 we will see an optimum path which has 
the control on the boundary for part of the path; the submarine travels 
at maximum depth for a while, then comes up. 

The subroutine for calculating u, v when w = w is contained in 

max, 

Appendix I, part C. The above features can be seen by examining either 
the flow chart for subroutine BOUNDW or the subroutine itself which is 
given in part G of Appendix I. 

Modifications similar to those made for w on the boundary would be 

made if v = v for some part of the route. One additional change re- 
max ^ ° 

quired in this case, ^ich was not necessary when w = w , is that ^ 

^ max 

be set equal to zero in the numerical routine whenever v = v and 
^ max 

< 0. V = and > 0 imply that the velocity is decreasing or 

moving away from the bound and hence the numerical routine described in 

equations (5.22) for generating corrections to the control variables 

would be used whenever this is the case. It should be noted that setting 

(Tw equal to zero in the case of'^w = w was not required since cTw does 
^ max 

not enter into the equations for the numerical routine in section 5. 

When V = V and H < 0, we get from equations 

max V 
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(6.4) 



H du + H dw = - H 
uu uw u 

H du + H dw = - H 

WU WW W 



the corrections to the control variables u, w as 



(6.5) 





- H 

u 


H 

uw 


du = 


- «w 


«ww 




\u 


«uw 




«wu 


«ww 




H 

uu 


- H 

u 


dw = 


«wu 


- H 

w 




Huu 


«uw 




«wu 


H 

WW 



Using the du, dw found in equations (6.5), we get equations for correct- 
ing u, w 



( 6 . 6 ) 



u„ = u , + du 
n n-1 n 



w_ = w , + dw 
n n-1 n 



If both V = V and w = w and it is also true that H 4 0 and 
max max v 

H < 0, we generate corrections to the control variable u by using the 
w 

equation 

(6.7) H du = - H 
uu u 



from which 

( 6 . 8 ) 



- H 



du = 



u 



H 



uu 



and the iterations to correct u take the form 



(6.9) u„ = u - + du . 

n n-1 n 

Both of these situations were encountered in generating Path IV 
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which will be discussed in section 8. The results given there will indi- 



cate that the modifications needed when v = v 



or V = V 



and w = w 



max 



max 



max 



do not affect convergence of the Newton- Raphson iteration schemes in the 
numerical routine, because we will find that Path IV is admissible. 

The flow chart for the subroutine BOUNDV, which is used when v = 
Vmax> seen in part D of Appendix I. Part B of Appendix I con- 

tains the flow chart of subroutine VUW which takes care of the case where 

V = V and w = w . The result of setting <Tv = 0 when v = v and 

max max max 

< 0 can be seen by examining the flow chart of the numerical routine 
in part A of Appendix I. The deck listings of BOUNDV, VUW, and the nu- 
merical routine can all be found in part G of Appendix I. 

7. Paths. 

Our computations have established the existence of three extremals. 
The three paths all satisfy the Euler-Lagrange equations as outlined in 
section 4. 

These three paths can best be compared by listing the contrasting 
points of the three. The areas in which the greatest difference appeared 
among the three routes are the probability of detection p(T) , the depth 
w, and the constants (h^,h 2 ) which yield admissibility. 

The following is a list of the results for each path in the three 
areas just mentioned. 



Path I 



p(T) 



15741 X 10 



-3 



w 



thermocline 200 feet 




(-.00019, .00024) 



Path II 



p(T) 



22725 X 10 



-3 



29 



w w up to 683 feet 

max 

(h^,h2) (-.00038, .00048) 

Path III 

p(T) .23060 X 10"^ 

w 525 to 725 feet 

(h^,h2) (-.00037, .00046) 

Path I gives us an absolute minimum, i.e. , a minimum under either 
weak or strong variations. Path II gives us a relative minimum if only 
weak variations are allowed. Path III is an extremal, but does not fur- 
nish a relative minimum under either weak or strong variations. We 
called Path III a worsimax path. 

If any additional information concerning any one of the above paths 
is desired, copies of the three paths can be found in Appendix II. Given 
there is a printout-, of each path with the coordinates (x,y) denoting the 
submariners location, the control variables u,v,w, and the probability 
of detection p(t) for each time step. 

Note that in Path II w = w for nearly all of the route. In cal- 

max 

culating this path, the iteration scheme described in equations (6.1), 
(6.2), (6.3) was used. The fact that Path II converges to the desired 
terminal point (x,^,,y,^) substantiates the results given in section 6. 

Analysis of Paths . Checks on the conditions as outlined in section 
4 point out the following facts. Path I satisfied the following condi- 
tions: 

1. admissibility conditions, 

2. Euler equations, 

3. Legendre conditions, and 

4. Weierstrass condition. 
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Path II satisfied 



1. admissibility conditions, 

2- Euler equations, and 
3. Legendre conditions, 
but not the Weierstrass condition. 

Path III satisfied 

1. admissibility conditions, and 

2. Euler equations, 

but not the Legendre conditions nor the Weierstrass condition. 

Hence, by the criterion established in section 4, there is but one 
minimum path, that being Path I. Note that the check of the envelope 
condition was not included in this investigation. 

The above results point out an important fact which is often ignored 
or overlooked; the generation of an extremal, by no means guarantees that 
you have the desired minimim. This emphasizes the need for a check on 
all of the conditions for a minimum at each time step of the numerical 
integration scheme for generating the path. If checking after each time 
step requires excessive computer time, the checks may be performed at 
some appropriate periodic intervals. 

It can be noted at this time, that if we restrict ourselves to weak 
variations as defined in section 4, both Path I and Path II are extremals 
which yield relative minima. In contrast to this, analyzing the paths 
and considering strong variations yields the result noted above, namely, 
that Path I is really the only relative minimum among the three paths. 

After Path II, which does not give a relative minimum, was genera- 
ted, a search, as outlined in the discussion of the Weierstrass condition 
in section 4, was used to determine a new set of control variables u,v,w, 
which would S£^isfy the Weierstrass condition. The resulting paths 
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turned out to converge to Path I. The subroutine search was used for 
this purpose and is contained in Appendix I. 

A similar search over the grid outlined in section 4 was performed 
when the Legendre conditions were not met in the case of Path III. The 
set of control variables which gave the smallest value for the Hamilton- 
ian in this search were then used to continue the numerical routine. 

Again the sequence of paths produced by the numerical integration con- 
verged to Path I. 

A path was judged admissible if it came within one-fourth mile of 
the desired endpoint; it was felt that further accuracy was not worth the 
computing time. 

To test the convergence of the numerical routine, on a few paths the 
routine was allowed to continue until no further improvement occurred. 

In each of the three paths above, duplication occurred with accuracy of 
at least one-tenth of a nautical mile. Duplication here means the abil- 
ity of the numerical routine to repeat itself after once converging to 
the desired terminal point. 

It has been noted that a spiral pattern of convergence about the 
desired terminal point is present in the computation of each path. It 
is not clear why this occurs but the following is offered as a possible 
explanation. In the computation, we take = 0 and vary u,v, 

w,x,y, h^jh^, but we assume that x,y have the values they assume on the 
path. For the purposes of explanation, let us examine the equation 
= 0 , from which we get an equation of the form 






H <^x + H <^y = 0 
ux uy ^ 
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when we take its variations. The last two terms in this expression drop 
out in linear problems and can be shown to be negligible if T is small 
in any case. They introduce considerable complication and extra computa- 
tion and hence were discarded. This omission may be the reason for the 
spiral pattern of convergence. It should be pointed out that the terms 
were omitted in only the correction routine. If the sequence of paths 
converges, there is no related error in the path to which they con- 
verge. 

In generating Path II, a subroutine was used within which the head- 
ing was fixed and a search over the depth and the speed was conducted to 
determine which set of values for these two controls gave the smallest 
value for the Hamiltonian'^H. These values were then used in the numeri- 
cal routine to insure a start in the proper direction. This method was 
put into use when it was found that poor initial choices for u,v,w caused 
the Newton- Raphson routine to diverge at the beginning of the route. This 
subroutine SEARCH can be seen in Appendix I, part G/ 

To insure a proper start in computing Path II, the subroutine WORSI 
was used. This subroutine is the same as subroutine BOUNDW described in 
section 6 with one exception, that being that w = is replaced by 

w = 500 or some other intermediate value for the depth. This subroutine 
then fixes w at 500 and computes u and v using the iteration scheme des- 
cribed in equations (6.3). The resulting set of values for u and v are 
then combined with w = 500 to make up the initial guesses for the numer- 
ical routine. This subroutine is also given in Appendix I, part G. 

Few problems were encountered in the generation of Path I, but the 

introduction of conditions to cause the submarine to assume its maximum 

depth w or velocity v aggravated the situation and introduced dif- 
^ max max 
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ficulties to make the subroutines listed above necessary. 

It should be realized that this problem contains some real difficul- 
ties, if approached blindly. With the proper background and forethought, 
most of the difficulties can be anticipated and handled, when encountered, 
by methods such as those described above. 

8. Corner. 

This section contains a discussion of a case in which the path gen- 
erated contains a corner. 

A corner appears when control variables u,v,w which minimize the 
Hamiltonian are discontinuous functions of t. For convenience, the nota- 
tion 




will be used to denote a set u,v,w of control variables. The conditions 

for a corner are, first, that there exists some point on the route, call 

-^1 -->2 

it t^, where two sets of control variables U and U both minimize the 

Hamiltonian, H, and second that one set of control variables, say, U , 

gives a smaller value for H when t ^ t^, whereas the other set of con- 
-^2 

trols U yield a smaller H for t> t^. These two conditions can be 
stated as, first 

(8.1) H(U^) = H(U^) = min t = t, 

1 



and second 

H(U^) < H(U^) for t ^ t 
( 8 . 2 ) 1 
H(U ) > H(U ) for t > 

for some neighborhood of t^. 

At a corner the numerical routine for the corrections should be 
amended by adding terms to handle the changes due to the variation of 
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the corner time t 



r 

The changes necessary are outlined in the earlier report -Optimum 
Submarine Routing’’ section 6. Even without these correction terms, 

the numerical routine generated an admissible path which contained the 
corner, namely Path IV* 

To construct a mathematical model in which a corner could be ex- 
pected the functions described in section 2 were changed in the following 
manner. We replaced 



1 




in the equation for f^ by 



with w^ equal to five hundred feet. The constants in f^, the equa- 

tions representing the passive and active defenses respectively, were 
chosen in such a way that a corner could be anticipated. 

The model used was one in which the passive defense was dominant 
for the beginning of the route, the two would become equal at approxi- 
mately the middle of the route, and the active search would then dominate 
for the latter part of the route. These facts are apparent when one looks 
at Path TV in Appendix II. 

When we analyze Path IV, we see that the submarine proceeds at ap- 
proximately thermocline depth for the first part of the route and then 

changes to w = w for the remainder of the route. It can also be noted 
® max 

that the speed, v, was considerably less than v = v until the corner 

max 
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was encountered, at which time v became equal to v . Whenever the active 

U13.X 

search completely dominates, the submarine goes as deep and as fast as 
possible. 

A check on the Hamiltonian, H, after each time step shows that H is 
constant, within the accuracy of the routine, for the time steps before 
we reach the corner, but is not constant as we proceed beyond the corner. 

It is not clear why H does not remain constant throughout the routej but 
a possible explanation is that the corner was effectively passed before 
it was found, i.e., the numerical routine failed to detect the corner 
when the conditions for a corner were in fact present. The failure to 
find the corner immediately is a result of the grid size used in the sub- 
routine SEARCH to compare the Hamiltonian for controls U generated by our 
numerical routine with the controls used in the search. As noted be- 
fore, this grid size was decided upon when a finer grid was found to re- 
quire excessive computer time. Considering that it takes a while for the 
routine to find the corner and it takes the routine a certain amount of 
time to settle down after the corner is found, the fact that admissibil- 
ity was accomplished was thought sufficient to justify omission of the 
corner correction terms. 

Notice that in Path IV v = v with w<w for a part of the route 

max max 

before the corner *>and then v = v and w = w for the portion of the 

max max ^ 

route that comes after the corner. Path IV was generated using the iter- 
ation schemes described in equations (6.4), (6.5), (6.6) when v = v 

max 

and w < w and by using (6.7), (6.8), (6.9) when v = v and w = w 

max max max 

The admissibility of Path IV confirms the results stated in section 6 for 
V on the boundary and when both v and w are on the bound of their allow- 
able values. 
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9. Observations. 



This section contains observations which may be helpful to a person 
wishing to continue the study of the submarine routing problem. 

In the Newton-Raphson iterations which occur it may be possible to 
improve both convergence and accuracy as follows. For example, in gen- 
erating the control variables let us make up an error function 



If each successive iteration does not diminish this error function, then 
we should diminish the preceeding corrections by a factor of say, two, 
or five. The reason is that the Newton-Raphson iteration moves the var- 
iables in the right direction but may overshoot if the linear terms are 
not dominant. The iteration would be terminated when the above error 
function was less than some preassigned value. The incorporation of such 
routines might well improve convergence, save computing time, and improve 
accuracy. The convergence criterion above is derived from the fact that 
satisfaction of the Euler equations implies H ,H , H are all equal to 
zero. Similar conditions could be established for the other Newton- 
Raphson iterations for generating corrections to the parameters h^,h 2 * 

In the iterations to correct these the established criterion would be 
based upon admissibility of the path. By letting (x,y) represent the 
endpoint of the path generated and be the desired terminal point, 

we could write the condition as 



2 



2 : 
+ e. H + e.’ll 
2 V 3 w 



2 




(x - x^)^ + (y - < £ 
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10. Summary. 

In this paper the submarine routing problem was studied. Functions 
were chosen that seemed to be typical of the functions representing the 
detection devices, both passive and active. Information that could be 
arrived at only through the use of empirical data such as the sea state, 
for example, was not made a part of these functions. 

Using this mathematical model and determining the path from (xQ,yQ) 
to (xrjxyYrjd 9 in a fixed time T, which minimizes the probability of detec- 
tion, p(T) , resulted in the generation of three extremals. Examination 
of the paths using established criterions for a relative minimum lead to 
the following results: one path yielded the desired minimum, a second 
satisfied all conditions except the Wfeierstrass condition and the third 
path was just an extremal, satisfying neither the Legendre nor the Weier- 
strass conditions. 

Situations were encountered in which the speed, v, or the depth, w, 
or both V and w were on the boundary of allowable controls. It was found 
that if the control on the bound was set equal to the boundary value and 
corrections generated for the remaining control variables, admissibility 
was accomplished just as when all controls were interior to their allow- 
able ranges. 

With a change in the original model for the active defense, condi- 
tions conducive to a corner were established. The numerical routine 
then generated a path which had a corner and was admissible. Admissi- 
bility here was accomplished without the use of corrections for the 
corner, and for this reason the corrections were not made a part of the 
numerical routine. 

The numerical routine as described in section 5 was programmed in 
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Fortran 1960 for a GDC 1604 computer. Both a flowchart and the program 
for this routine can be found in Appendix I, the flow chart in part A 
and the deck listing for the program in part G. 
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Appendix I 



This appendiic contains the f]^ovj charts and the 
deck listings for the program of the numerical routine 
and the subroutines that were used to generate the 
paths described in the text. 

A, Flow chart for the numerical routine. 



.^k 



Dimension mRS,YC,DI,C,AK,TAU,X,Y,Z,UT,VT,WT 

t 

A 



/^ad XT , T . THETA , A1 , B1 , Cl , Dl , SI , C2 , 
D2 ,X2 , Y2 ,H1 ,H2 ,FAC ,FMUL,V/0 ,LMAX , 






'Q = THETA / 57.29 

COS = COSF(Q) 
SIK = SIl^(Q) 
,X2 = SI COS + XT 
I Y2 = SI SIN 
; VAV = XT / T 



...Mi- 



Print XT ,T ,TEETA ,A1 ,B1 , Cl ,D1 ,S1 . C2 
D2 ,X2 , Y2 , HI e H2 ,FAC ,FMOL,WO ,LM4X , 

V wl .EMiX.VIiyL 



I C(l) = 0.0 
! C(2) = 0.0 
, C(3) = 0.0 
. C(k) = 0.0 
ixSTSP = T / FAC 



N1 = f / XSTEP i 
XNl = N1 
STEP = T / XNl 
N2 = N1 + 1 



i 

\ 

\/ 



!u = - .21 
, V = VAV I 
' V/ = wo 1 



41 



0 



rri^n 

1 ' 





TYPJR 


'4 

= 0.0 







lu = UZERO 
= VZERO 
|W = V/ZERO 




N 

YVARSCl) = 


0.0,1 = 1,9 1 


S 


/ 




Compute DETER,mil,FU>12>F^ 
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UZiiRO" = if' 
VZERO = V 
WZERO = ¥ 



rCompute FIl ,F I2>r , FY,Sll. S12 , S2_l , S2^- 
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( Print FROB) 



STOPl 
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3, Flow chart for subroutine VUV/. 
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C» Flow chart for subroutine BOUtJDV/, 



START! 

v!/... 

W = 



\,/ 

ICompute HU.HV.KUU.HVU.HW 

- HVU HVU ' 
GUI41 = - HU HW + HV HVU 
GUI^I 2 = - Hy HUU HU HVU 








|dU = GUlIi r DEN 
|DV = GUl'12- / DEN 
1 

^ 

U = U + DU 

V = V + DV 



I Compute HU ,HV,HUU, HVU.HWl 

i 

i 

IDSN = EUU HW - HVU HVUl 
jGlDfl. = - HU HW + HV HVU 
iGIMB = .. HV HUU HU HVUl 
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D 



'loiv chart for subroutine BOUNDV 



I START 






V = 


VMA.XI 


V 


/ 



VI 



Compute HU, HV J.HDU,ffi^.HWW 



IDEN = HUU Hlfw' - HItIU HWU 
Groil = - HU H*JW + HW EVRJ 



GU> 1 2 = - m -J HUU + HU ErfU 
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Flow chart for subroutine SEARCH, 



i START I 






! Dimension VIKT,VflCNT,H I 



-V- 



HAM = XLAl-I V COSU + Xl^ V SINU + FIl'+ FI2| 



V/ 




t ; 




E^Il 



\l/ 



N = 1 
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F 



Flov; chart for subroutine WORSI 




iu = U + DU 
!v = V + DV 




\ 


/ 




Compute HU,HV,HUU,HVU,H\AM 



^ 

DEN = HUU HW - HVU H7U 
GUMl = - HU HW + HV HVU 
GUM2 = - HV HUU + HU HVU 




Go. Fortran Program: Printout of deck of punch cards 

submitted to computer. 



The main nimerical routine and various subroutines 
are given here. The equations in the ntimerical routine 
and the first fovee subroutines as listed here are the 
ones which gave a corner. The equations in subroutine 
k’ORSI are the ones used in generating the three extremals. 



PROGRAM SUBROUTE 

C YVARS ( i )=LAM3 YVARS(2)=MU3 YVARS(3)=A11 YVARS(4)=A12 YVARS(5)*A21 

C YVARS(6)=A22 YVARS(7)= Z YVARS(8)= X YVARS(9)= Y 

C XLAM=H1+LAM3 XMU=H2+MU3 

DIMENSION YVARSC9) » YC ( 9 ) » DY ( 9 ) ,C ( 4 ) » AK ( 4 » 9 ) »TAU(400) »X(900) ♦ 

+ Y(90 ), 2(400), VT(400) »UT(400) ♦WT(400) 

READ 1 »XT,T,THETA,A1 ,B1 ,C1 ,D1,S1 »C2,02,X2,Y2 ,H1 »H2,FAC»FMUL, 

+WO , UMAX » W 1 , WMAX , VMAX 

1 FORMAT (-5.0,2F4.0,F2.0,F3.2,F2.1»F5.4,F4.0,F7.6,F5.4»F6.1»F4,1, 

+ F6«4,F6.0,F4.0,F4.2,F4.0, I2/F4.0,F5.0,F3.0) 

Q=THETA/57.2957 7951 
C0S=C0SF(0) 

SIN=SINF(Q) 

X2=S1*C0S+XT 
Y2 = Sl-«-SIN 
VAV=XT/T 
PRINT 2 

2 FORMAT ( 1H02HXT4X1HT3X5HTHETA2X2HA12X2HB12X2HC12X2HO12X2HS12X4HWMAX 
+1X4HVMAX) 

PRINT 3,XT,T,THETA,A1,B1,C1,D1 ,S1 , WMAX , VMAX 

3 FORMAT (3F5.0,F5,1»F4.2»F4.1,F5.4,F5.0,F5.0,F3.0/) 

PRINT 4 

4 FORMAT! 3H C24X2HD24X2HX24X2HY25X2HH17X2HH23X3HFAC2X4HLMAX 
+3X4HFMUL2X2HW05X2HW1 ) 

PRINT 5»C2»D2,X2,Y2,H1 ,H2 , FAC , UMAX ,FMUL ,W0,W1 

5 FORMAT ( F6 • 5 , F5 • 1 , F6 • 1 , F4 ♦ 1 , 2 F9 • 5 , F6 • 0 » 1 3 » F8 • 2 , F6 *0 , F6 • 0/ ) 

C DEFINE RUNGE-KUTTA CONSTANTS. 

C(1)=0.0 



I 

I. 
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C(2)=0.5 
C(3)=0.5 
C(4)=1.0 
X(1)=0.0 
Y( 1 )=0.0 
TAU(1)=0.0 

C COMPUTE TIME STEP LENGTH. 

XSTEP=T/FAC 

N1=T/XSTEP 

XN1=N1 

STEP=T/XN1 

N2=N1+1 

C GUESS INITIAL VALUES FOR U, V, AND W BEFORE FIRST ITERATION 
U=-.9 
V = VAV 

w=wo 

C ENTER LOOP FOR GENERATING CORRECTIONS TO H1»H2. 

DO lA L=1,LMAX 
TVAR=0.0 

C GUESS INITIAL VALUES FOR U»V»W» AFTER FIRST ITERATION 
IF (L-1) 198»199»198 

198 U=UZERO 
V=V2ERO 
W=WZERO 

199 UT(1)=U 
VT ( 1 )=V 
WT { 1 ) =W 

DO 6 1=1,9 
5 YVARS( I )=0,0 

C enter LOOP FOR COMPUTING THE PATH FOR EACH TIME STEP. 

DO 11 K=2,N2 
IF(K-3) 901,902,902 

902 CALL SEARCH! A1,31,C1 ,D1 ,FI1 ,FI2»XLAM,XMU,C0SU,U,V,W,SINU,W0, 
+G1,G2,W1,C0SUMQ) 

C ENTER LOOP FOR RUNGE-KUTTA INTEGRATION. 

901 DO 88 1=1,4 
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DO 7 J=l,9 

7 YC( J)=YVARS( J)+C( I )*AK( 

XLAM=H1+YC( 1 ) 

XMU=H2+YC(2 ) 

XD=(XT-YC(8) )*C0S+S1-YC(9)*SIN 
FI1=EXPF(-(XD*L06F(2.)/1000.) ) 

FI2=C2+D2*EXPF(-( (YC(8 ) -X2 ) **2+ ( YC ( 9 ) -Y2 ) **2 ) *LOGF ( 2 . ) / 100000 . ) . 
COSU=COSF(U ) 

5liNU=SINF(U) 

C0SUiMQ = C05F(U-Q) 

SINUMQ=SINF(U-Q) 

G i = C 1 + D 1* ( W-WO ) * ( W-WO ) 

G2 = 1 .+0 1* ( W-WO ) * ( W-WO ) 

63= ( .01 )*FI 1*(A1+B1*V*V) 

6A=(-61*(2o*Dl*(W-W0) )+62*(2.*Dl*(W-W0) ) )/(G2*62) 

65=(Cl+Dl->:-( W-W 0 ) * ( W-WO ) ) / ( 1 . +D 1 * ( W -WO ) * ( W-WO ) ) 

67 = 1.+W*W/(W1---W1) 

Pl = ( ( i .+01--=- (W-WO) W-WO) )*2.-«-Dl+4.*Dl»( W-WO) *D1*( W-WO) ) / ( G2*G2 ) 
P2=(-( 8.->Dl*(W-W0)*( W-WO) )*D1) /(G2^:-G2 ) 

P3= (-( Cl+Dl^:-( W-WO )*( W-WO) )*2.*D1-4.*D1* (W-WO) *D1*( W-WO) ) / (62*G2) 
P4= (G1«-8.*D1-K-(W-W0)*D1*(W-W0) ) /(62*G2*G2) 

P9 = FI 2* (-2 . / ( W1*W1^^-G7*G7 ) +8 .*W*W/ ( W1*W1^«-W1^^W1*G7«-G7*67 ) ) 
HV=XLAM*COSU+XMU*SINU+ ( .01 ) *F 1 1'-- ( 2 . *B1*V ) *G5* ( 1 .-.2 5*COSUMQ ) 
HU=-XL AM*V*S I NU+XMU*V*C0SU+G3^^65* ( . 2 5*S I NUMQ ) 

HVV=( .01 )^^2.*B1--^FI l->G5*(l.-.25«-C0SUMQ) 

HVU=-XLAM*SlNU-rXMU*COSU+G5*.Ol*FIl*2.*Bl*V*{ .25*5 I NUMQ) 
HUU=-XLAM*V*C0SU-XMU*V*S I NU+G3*G5* ( . 25*COSUMQ ) 

HW = G3*( l.-.25*(C0SUMQ) ) *G4+F I 2* ( -2 .* ( W/Wl ) * ( 1 . /W1 ) ) /(G7*67 ) 
h'WU=63*( .25*(SINUMQ) )*64 

HWV=G4*( .01 )*FI 1*( 2.*B1*V)*( l.-.25*(COSUMQ) ) 

KWW= 63*( ( l.-.25*COSUMQ)*(Pl+P2+P3+P4) )+P9 

P5=HVV*HWW-HWV*HWV 

P6= (-HV)*HWW+HW*HWV 

P7=HVV*(-HW )+HWV*HV 

DETER=HUU*P5-HVU* ( HVU*HWW-HWU*HWV )+HWU* ( HVU*HWV-HWU*HVV ) 

FUM1= ( -HU ) *P5-HVU*P6+HWU* ( ( -HV ) *HWV+HW*HVV ) 
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c 

c 



FUM2=HUU^^-P6+HU---( HVU^-HWW-IHWU^HWV ) +HWU'’'^ ( HVU* ( -HW ) +HWU*HV ) 

FUM3 = HUU-*i-p7-H\/U*( HVU >«• ( -HW ) +HVJU*-HV ) -HU* ( HVU*HWV-HWU*HVV ) 
rAi I furroUTINE FOR CORRECTIONS TO U*V»W» 

CALL VUW( deter ♦ FUMl » FUM2 » FUM3 ♦ XLAM*XM'J ,FI1»FI2»A1»B1*C1*D1»W0» 



+W1 ,WMAX»VMAX»U»V»W»Q) 

FIV=COEFF OF FIl F2V=COEFF OF FI2 

IF(K-2) 195»297»19.5 

297 IF (I-l) 195»196*195 
196 UZERO=U 
VZERO=V 
WZERO=W 

195 XD=(XT-YC( 8 ) ) *COS+Sl-YC ( 9 ) *S I N 
FI l = EXPF(-(XD*LOGF(2;. ) /l000. ) ) 

FI 2=C2+D2*EXPF ( -{ ( YC ( 8 )-X2 )**2+ ( YC ( 9 ) -Y2 ) **2 ) *LOGF ( 2* )/ 100000. ) 
G1=C1+D1*(W-W0)*( W-WO) 



G2 = 1 . +D 1* ( W-WO ) * ( W-WO ) 

G3=( .01)*FI1*(A1+B1*V*V) 

GA= (-Gl*( 2.*D1*(W-W0) )+G2*(2.*Dl*{ W-WO) ) )/(G2*G2) 
G5=(C1+D1*(W-W0)*(W-W0) )/( l.+Dl*(W-WO)*(W-WO) ) 



G7=1.+(W/W1)**2 

COSU=COSF(U) 

SINU=SINF(U) 

SINUMQ=SINF (U-Q ) 

Pl=^a^Sl* (V/-WO)*( W-WO) )*2.*D1+4.*D1*( W-WO) *D1*( W-WO) )/(G2*G2) 

P2= { -( 8.*D1*( W-WO) *( W-WO) ) *D1 ) / ( G2*G2 ) 

P3= ( - ( Cl+01* ( W-WO ) * ( VJ-V/0 ) )*2,*D1-4.*01*( W-WO ) *D1* ( W-WO ) ) / ( G2*G2 
pz).= ■ ( G1*8 .*D1* ( W-WO ) *D1* ( V/-WO ) ) / ( G2*G2*G2 ) 

P9=FI 2* ( -2 • / ( W1*W1*G7*G7 )+8.*W*W/ ( W1*W1*W1*W 1*G7*G7*G7 ) ) 
HV=XLAM*COSU+XMU*S INU+ ( .01 ) *Fl 1* ( 2 .*B1*V ) *G5*( 1 .-.25*COSUMQ ) 
HU=-XLAM*V*SINU+XMU*V*C0SU+G3*G5*( .25*SINUMQ ) 



HVV=( ,01)*2.*B1*FI 1*G5*(1 .-.25*C0SUMQ) 

HVU=-XLAM*SINU+XMU*COSU+G5*.Oi*FI 1*2.*B1*V*( .25*SINUMQ ) 
HU'J=-XLAM*V*COSU-XMU*\/*S I NU+G3*G5* ( . 25*COSUMQ ) 

HW=G3*( 1.-. 25* ( COSUMQ ) )*G4+FI2* (“2 .*(W/W1) *( l./Wl ) ) /(G7*G7) 
HWU=G3*( .25*(SINUMQ) )*G4 



.1 



HWV=G4^-( .01 2.*31*V)*( l.-.25*(GOSUMO) ) 

HWW= G3*( ( l.-.25*COSUMQ)* (P1+P2+P3+P4) )+P9 
F1V=( A1+S1*V*V)*( .01 )-^ (l.-.25*COSUMQ)*Gl/G2 
F2V = l./(l.+ (W*W/(Wi^Wl ) ) ) 

FY = F1V*FI l*SIN-»LOGF( 2. ) / 1 000 . + F2V* ( -D2*EXPF ( - ( (YC (8 ) -X2)**2 + 

+ (YC(9)-Y2)**2)*LOGF( 2. )/100000. ) * ( 2 . * ( YC ( 9 ) -Y2 ) *L0GF ( 2 . ) /1 00000. 



+ ) ) 

FX = F1V*FI l*SIN-x-L0GF( 2. ) / 1000 . + F2V* ( -D2*EXPF ( - ( (YC(8)-X2)**2 + 

+ CYC(9)-Y2)-^*2)*LOGF( 2. )/100000. ) * ( 2 . * ( YC C 8 ) ~X2 ) *LOGF ( 2 . ) / 1 00000. 
+ ) ) 

DENOM=HVU* ( HVU*HWW-HWU*HWV ) -HUU* ( H VV*HWW-HWV*HWV ) +HWU* { HVV*HWU- 
+HWV*HVU ) 

511= ( ( HVV*HWW-HWV*HWV) «-(-V*SINU) + ( HVU*HWW-HWV*HWU )*( -COSU ) ) /DENOM 
Si 2= ( ( HVV*HWW-HWV*HWV ) * ( V^COSU ) + ( H VU*HWW-HWV*HWU ) * ( -S I NU ) ) / DENOM 
521 =( ( HVU^HWW-HWU*HWV ) * ( V*S INU ) + ( HUU*HWW-HWU*HWU ) »COSU ) /DENOM 
522 = ( ( HVU«-HWW-HWU*HWV ) * ( ( -V ) *COSU ) + ( HUU*HWW-HWU*HWU ) *5 1 NU ) / DENOM 
DY(1)=-FX 
DY ( 2 )=-FY 



I F ( V-VMAX )118,119»119 
118 DY ( 3 )=C0SU*S21-V*SINU*S11 
DY(4)=COSU*S22-V*SINU*S12 
DY( 5 )=SINU«-S21+V*C0SU*S11 






120 



DY ( 6 ) =5 INU*S22+V*COSU*S12 
GO TO 120 
V=VMAX 
DY ( 3 ) = 

DY( 4)= 

DY ( 5 ) = 

DY(6)= 



-V^-SINU*S11 

-V*SINU*S12 

V«-C0SU*S11 

V*C0SU*S12 



DY ( 7)=FI i*FlV+FI2*F2V 
DY(8 )=V*C0SU 



DY(9)=V*SINU 



DO 8 J=l»9 
3 AK( I »J)=STEP*DY(J) 
88 CONTINUE 
DO 9 J=l»9 
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9 YVARS( J)=YVARS( J)+(AK( 1,J)+2.*AK(2»J)+2.*AK( 3, J)+AK(4, J) )/6. 
TVAR=TVAR+STEP 
TAU(K)=TVAR 
X(K)=YVARS( 8) 

Y( K ) =YVARS( 9 ) 

Z(K)=YVARS(7) 

UT(K)=U 
VT(K)=V 
WT( K)=W 

11 CONTINUE 

DET=YVARS ( 3 ) *YVARS ( 6 ) -YVARS ( 5 ) *YVARS ( 4 ) 

XNUM1= ( XT-YVARS ( 8 ) ) *YVARS ( 6 ) +YVARS ( 9 ) *YVARS ( 4 ) 

XNUM2=-YVARS ( 9 ) *YVARS ( 3 ) - (XT-YVARS ( 8 ) ) *YVARS ( 5 ) 

H1=H1+ FMUL*XNUM1/DET 

H2 = H2+ FMUL*XNUiM2/DET 

PRINT 12 

12 FORMAT (2X2HN25X2HH17X2HH25X5HX(N2)4X5HY(N2)6X5HZ(N2) ) 

PRINT 13»N2»H1,H2»X( N2 ) »Y( N2 ) ,Z (N2 ) 

13 FORMAT (I4,2F9.5»2F9.1,E13.5/) 

18 PRINT 19 

19 FORMAT ( 1H03X3HTAU6X1HX7X1HY7X1HU5X1HV5X1HW10X1HZ/) 

PRINT 20, ( TAU( I ) ♦ X (I ) , Y ( I ) ,UT ( I ) , VT ( I ) ♦ WT ( I ) , Z ( I ) ♦ I » 1 ,N2 ) 

20 FORMAT ( F8 . 2 , 2 F8 . 1 , 3 F6 . 1 , El 3 . 5 ) 

■ 14 CONTINUE 

PROB = 1.0-EXPF(-0.0001*Z(N2 ) ) 

PRINT 21,PROB 

21 FORMAT (1H02X5HPROB=E13.5/) 

STOP 

END 



SUBROUTINE VUW ( DET ER , FUMl , FUM2 , FUM3 , XLAM ,XMU , F 1 1 , F 1 2 , A1 , Bl ,C1 ,D1 , 
+WO , W 1 , W MAX , VMAX , U , V » W , 0 ) 

DO 5 1=1,7 
DU=FUM1 /DETER 



DV=FUM2/DETER 
DW=FUM3/DETER 
U = U-rDU 
V=V+DV 
W=W+DW 

IF(WMAX -W) 390,390,392 
390 IF(VMAX-V) 1000,1000,1001 
392 IF(VMAX-V) 3000,3000,3001 

. 3000 CALL 30UNDV(U,V,W,XLAM,XMU,FI1,A1,B1,C1,D1,W0,FI2»Q,VMAX,HV,W1) 
IF(HV-O.O) 4000,4000,3001 
4000 RETURN 

1000 DO 2000 1=1,5 
V=VMAX 
W=V.'MAX 

COSU=COSF(U) 

SINU=SINF(U) 

C0SUKQ=C0SF(U-0) 

SINUMQ=SINF{U-Q) 

61=C1+D1*(W-W0)*{ W-WO) 

62 = 1 « +D 1* ( W-WO ) * { W-WO ) 

63=(..01)-^FI l^:-( A1+B1*V->^V) 

G4=(’-Gl*(2.^Dl->«-(W-WO) )+62*(2.*Dl*{ W-WO) ) )/(G2*G2) 

G5= (C1+D1*{ W-WO )*( W-WO) )/ ( 1,+D1*{ W-WO)*- (W-WO) ) 

HU = -XLAM*V*SINU+XMU*V*C0SU+G3*-G5*-( .25*SINUMQ) 

HUU = -XLA.M*V*COSU-XMU*V*SlNU+63*G5-5f { *25*COSUMQ) 

DU=-HU/HUU 

U=U+DU 

2000 CONTINUE 
RETURN 

1001 CALL BOUNDW(U,V,W,XLAM,XMU,FIl ,A1,B1,C1,D1,W0,FI2,W1,Q,WMAX,HW) 
IF(KW-O.O) 5000,5000,3001 

5000 RETURN 
3001 COSU=COSF(U) 

SIHU=SINF(U) 

COSUMQ=COSF(U-Q) 

SINUMQ=SINF(U-Q) 



1 

I 

) 
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G1 = C1+D1*(W-W0) *{ W-WO). 

G2 = l.+Dl*(W-V/0)*( W-WO) 

G3=( .01 1* ( A1+31*V^;-V) 

G4= (-Gl*( (W-WO) ) +G2* (2.*Dl* ( W-WO ) ) ) / (G2*G2 ) 

G5=(C1+D1*(W-W0)*( W-WO) )/( l.+Dl*( W-WO)* (W-WO) ) 

G7 = 1.+W*W/{W1«-W1) 

Pl=( ( l.+Dl*(W-WO)*(W-WO) )*2.*D1+4.*D1*(W-W0)*D1*(W-W0) )/(G2*G2) 

P2 = ( - ( 8.*D1* ( W-WO) *( W-V.'O) ) *D1 ) / (G2*G2 ) 

P3= ( - ( Cl+D 1* ( W-WO ) * ( W-WO ) ) *2 ,*D1-4.*D1* ( W-WO ) *D1* ( W-WO ) ) / { G2*G2 ) 
P4= ( Gl*8. *D1*( W-WO )*D1*( W-WO) ) /(G2*G2*G2) 

P9=FI2*(-2./ ( W1*W1*G7*G7) +8.*W*W/( W1*W1*W1*W1*G7*G7*G7) ) 
HV=XLAM*COSU+XMU*SINU+( .01 ) *F 1 1 * ( 2 . *B1*V ) *G5# ( 1 . 25*COSUMO ) 
HU=-XLAM*V*SINU+XMU*V*C0SU+G3*G3*( .25*SINUMQ ) 

HVV=( ,01)*2.*B1*FI 1*G5*( 1.-.25*C0SUM0) 

HVU=-XLAM*SlNU+XMU*COSU+G5*.Ol^tFI 1*2.*B1*V*{ .25*SINUM0) 
HUU=-XLAM*V*C0SU-XMU*V*SINU+G3*G5*( .25*C0SUM0) 

HW = G3*( l.-.2 5*(COSUMQ) ) *G4+F I 2* ( -2 .* ( W/Wl ) * ( l./Wl ) ) / (G7*G7 ) 
HWU=G3*( .25*(SINUMQ) )*G4 

HWV=G4*-( .01 )*FI 1*( 2.*B1*V)*( l.-.25*(COSUMQ) ) 

HWW= G3*( ( l.-.25*COSUMO)*(Pl+P2+P3+P4) )+P9 

P5=HVV*HWW-HWV*HWV 

P6= ( -HV ) *HWW+HW*-HWV 

P7=HVV* ( -HW ) +HWV*HV 

DETER=HUU*P5-HVU* ( HVU*HWW-HWU*HWV ) +HWU* ( HVU*HWV-HWU*HVV ) 

FUM1= ( -HU ) *P5-HVU*P6+HWU* ( ( -HV ) *HWV+HW*HVV ) 

FUM2=HUU-x-P6+HU^^ ( HVU*HWW-HWU*HWV ) +HWU* ( HVU* ( -HW ) +HWU*HV ) 
FUM3=HUU*P7-HVU*( HVU* ( -HW ) +HWU*HV ) -HU* ( HVU*HWV-HWU*HVV ) 

5 CONTINUE 
RETURN 
END 
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SUBROUTINE BOUNDW ( U » V , W » XLAM » XMU » F 1 1 , A1 , B 1 » C 1 » D 1 » WO » FI 2 » W1 »Q*WMAX» 
+HW) 

W=WMAX 

COSU = COSF(U ) 

SINU=SINF(U ) 

COSUMQ=COSF (U-Q) 

SINUMQ=SINF (U-Q) 

G1=C1+D1*(W-W0) *{ W-WO) 

G2 = l»+Dl*(W-WO) W-WO) 

G3=( *01)*FI l^i-{Al + Bl-"‘V*V) 

G4= (-Gl*( 2.*Dl*(W-WO) )+G2*(2.*Dl*(W-WO) ) )/(62*G2) 

G5=(C1+D1*{ W-WO)*( W-WO) )/ ( l.+Dl*( W-WO)*(W-WO) ) 

G7=1«+(W/W1 )**2 

HV = XLAM->C0SU+XMU*SINU+( .01 )*FI 1 *( 2.*Bl*V)*G5^^n.-.25*COSUMO) 

HU = -XL AM^-V*S I NU+XiMU^-^V^COSU+GS^^GS-^ { . 2 5*S I NUMQ ) 

HVV=( .01)*2.*B1*FI 1*G5*( 1 .-.25^^C0SUMQ) 

HVU=-XLAM*SINU+XMU*COSU+G5*.01*FI1*2.*B1*V*( .25*SINUMQ) 
HUU=-XLAM*V*C0SU-XMU*V*S I NU+G3-»G5* ( . 25*C0SUMQ ) 

HW=G3*( i.-.25«-{COSUMQ) ) *G4+F 1 2* ( “2 .* ( W/Wl ) * ( l./Wl) )/(G7*G7) 

DEN=HVU*HVU-HUU*HVV 

GUM1=-HV*HVU+HU*HVV 

GUM2=HVU*(-HU)+HV*HUU 

DO 981 1=1,7 

DU=GUM1/0EN 

DV=GUM2/DEN 

U=U+DU 

V=V+0V 

COSU=COSF(U ) 

5INU=SINF(U ) 

COSUMQ=COSF(U-Q) 

SINUMQ=SINF(U-Q) 

G1=C1+D1*( W-WO ) * ( W-WO ) 

G2 = l.-rDl*(W-WO) W-WO) 

G3=( .01)*FI1^-(A1+B1*V>-V) 

G4={-G1*{ 2.*D1*(W-W0) )+G2^(2.*Dl*( W-WO) ) )/(G2»G2) 

G5= { C1+D1*{ W-WO) «•( W-WO) )/( l.+Dl*( W-WO)* (W-WO) ) 



HV=XLAM*COSU+XMU*SINU+ ( .01 ) -^F I 1* ( 2 . *B1* V ) *G5* ( 1 . - . 2 5*COSUMQ ) 
HU=-XLAM*V*5INU+XMU«-V«-COSU+G3*g5*{ .25*SINUMQ ) 

HVV=( .01 )-)^2.*Bl"-FI 1*G5*( 1 .-.2 5«C0SUMQ) 

HVU = -XLAM*SlNU+XMU*COSU+G5*.Ol^^FI 1*2.*B1*V*( .25*SINUMO) 
HUU=-XLAM^ V<-COSU-XMU*V*S I NU+G3*G5* ( . 2 5*COSUMQ ) 
DEN=HVU*HVU-HUU^f-HVV 
GUM1=-HV*HVU+HU*HVV 
GUM2=HVU*(-HU)+HV*HUU 
981 CONTINUE 
RETURN 
END 



SUBROUTINE BOUNDV ( U ♦ V , WtXLAM » XMU , F 1 1 , A1 , B1 »C 1 . D1 . WO ♦ FI 2 ♦Q» VMAX , 
+Wi ) 

V=VMAX 

COSU=COSF(U ) 

SINU=SINF(U) 

COSUMQ=COSF(U-Q) 

SINUM0=SINF(U-Q) | 

G1=C1+D1*(W-W0)*(W-W0) 

G2=1.+D1*(W-W0)*{ W-WO) 

G3= ( .01 )*FI 1*(A1+B1*V*V) 

G4= (-Gl*( 2.^^D1* (W-WO) )+G2*(2.*D1*(W-W0) ) )/(G2*G2) 
G5={Cl+Dl*(W-WO)*{>y-WO) )/{ l.+Dl*(W-WO)*(W-WO) ) 

G7 = l.+W^fW/ ( W1*W1 ) 

Pl=( { l.+01*(W-W0)^'^(W-W0) )*2.*D1+4.*D1*(W-W0) *D1*{ W-WO) )/(G2*G2) 
P2 = ( - { 8 . *D 1 * ( V/'-WO ) ( W- WO ) ) *D1 ) ’/ { G2 *G2 ) 

P3= (-(C1+D1*(W-W0)*( W-WO) )^^2.^^D1-4.*D1*(W-W0)*D1*(W-W0) )/(G2*G: 
P4= (G1*3.->D1^:-(’W-W0)*D1*( W-WO) ) /(G2*G2*G2) 

P9 = F I 2* ( -2 . / { W l*Wl-i«-G7*G7 ) +8 .*W^«-W/ ( W1»W1*W1*W1*G7#G7*G7 ) ) 
HU=-XLAM*V^^SINU+XMU*V-"-C0SU+G3*G5*( .2 5*SINUMQ) 
KUU=-XLAM^^V*COSU-XMU«-V*S I NU+GS^GS^^- ( . 25*COSUMQ ) 

HW = G3<-( l.-.2 5--(COSUMQ) ) *G4+F I 2* ( “2 . * ( W/Wl ) * { l./Wl ) ) / (G7*G7 ) 

HWU = G3-- ( .25*(SINUMQ) )*G4 
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HWW= G3*( ( 1,-.2 5*C0SUMQ)^«-(P1+P2+P3+P4) )+P9 

DEN=HWU*hWU-HUU*HWW 

GUM1=-HW*HIVU+HU*HWW 

GUM2=-HU^-HWUtHW*HUU 

DO 981 I=l*7 

D0=GUM1/DEN 

DW=GUM2/DEN ■ ' 

U=U+DU 

W=W+DW 

COSU=COSF(U ) 

SIKU=SINF(U ) 

COSUMQ=COSF(U-Q) 

SINUMQ=SINF(U-Q) 

G1=C1+01*( W-WO) *( W-WO) 

G2 = l. + Dl*(W-WO) '»( W-WO) 

G3=( .01)*FI l-K-( A1+B1*V*V) 

G4=(-G1*(2.*D1*(W-W0) )+G2*(2.<-Dl*(W-WO) ) )/(G2*G2) 
G5=(C1+D1*(V/-W0)*( W-WO) )/ ( l,+Ol*( W-WO)*(W-WO) ) 

G7 = l.+W*W/( V/1*W1 ) 

Pl = ( ( 1 .+Dl-^(W-WO)*( W-WO) )*2.*D1+4.*D1*(W-W0)*D1*( W-WO) ) / { G2*G2) 

P 2 = ( - ( 8 . -x-D 1 * { W- WO ) *.{ W- WO ) ) *0 1) / { G2 *G2 ) 

P3= (-{ C1+D1^^-(W-W0)*( W-WO) )*2.*D1-4.*D1*(W-W0)*D1*(W-W0) )/(G2*G2) 
P4= (Gl*8o*Dl^^-( W-W0)*D1*(W-W0) ) /(G2*G2*G2) 

P9 = FI2*(-2, / (Wl-x-Wl*G7^:-G7)+8.-^^W*W/ {Wl-i^Wl*Wl*Wl^^G7»G7*G7) ) 

HU = -XLAM-x-V-xs I NU+XMU*V*C0SU+G3*G5* ( . 2 5*S I NUMQ ) 
HU’U=-XLAM-x-V^^-COSU-XMU-xV-x-SINU+G3'X-G5*{ .25*C0SUMQ) 

HV;=G3--( l.-.25*(COSUMQ) ) ^G4+F I 2* ( “2 ( W/Wl ) * ( 1 . /W1 ) ) /(G7*G7) 

HWU=G3*(.25*(SINUMQ) ) *G4 

HWW= G3"x-{ ( l*-.2 5*COSUMQ) *(P1 + P2+P3+P4) )+P9 
DEN=HWU*KWU-HUU*HWW 
GU.Ml=-HW*HWU+HU^fHWW 
GUK.2=-HU-x-HWU-rHW*HUU 
981 CONTINUE 
RETURN 
END 



t 
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SUBROUT INE SEARCH ( A1 , B 1 , Cl , D1 , F 1 1 » F I 2 » X LAM tXMU » COSU . U . V , W ♦ S I NU , 
+W0»G1»G2»VJ1 , COSUMO ) 

DIMENSION VINT(l0j,WINT(10),H(50,50) 

F1V=( A1+B1*V*V) *( ,01 ) * ( l.-.25*COSUMQ)*Gl/G2 
F2V=1./(1.+(W*W/(W1*W1 ) ) ) 

HAM =XLAM*V*C0SU+XMUW*.SINU + F1V*FI1+F2V*FI2 

VINT(0)=0,0 

WINT(0)=0.0 

DO 131 M=l,5 

VINT (M) =VINT ( M-1 ) +3.0 

DO 122 N = 1 , 5 

WINT(N) =WINT(N-l)+200.0 

G1=C1+D1*(WINT(N)-W0)*(WINT(N)-W0) 

G2=1.+D1*(WINT(N) -WO ) * ( W I NT ( N ) -WO ) 

F1V=(A1+S1->^VINT(M-)*VINT(M) )*( ,0l)*( 1 . - . 2 5*C0SUMQ ) *G 1 /G2 
F2V = 1./(1,+WINT(N)*WINT(N)/(W1-»W1) ) 

H(M,N) =XLAM*VINT( M)*COSU+XMU*VINT( M)*S INU+F1V*FI 1+F2V*FI2 

122 CONTINUE 
131 CONTINUE 

HAMSRCH=H( 1 ♦! ) 

VSRCH=VINT( 1) 

WSRCH=WINT( 1) 

DO 121 M=l»5 
DO 123 N=l,5 

IF (HAMSRCH-H(M,N) ) 123»123»lll 
111 HAMSRCH=H(M»N) 

VSRCH=VINT(M) 

WSRCH=WINT( N) 

123 CONTINUE 
121 CONTINUE 

IF (HAMSRCH-HAM) l»2i2' 

1 V=VSRCH 
W=WSRCH 

2 RETURN 
END 
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SUBROUTINE WORSI ( U » V » W ,XL AM » XMU ♦ F 1 1 » A1 . B1 » C 1 » D 1 , WO ♦ F 1 2 ♦ W1 *0 ) 
W=500. 

SINU=SINF(U) 

COSUMQ=COSF(U-Q) 

SINUMQ=SINF (U-Q) 

G 1 = C 1+D 1* ( W-WO ) -X- ( W-WO ) 

G2 = l. + Dl'X(W-W0)-K-( W-WO) 

G3= ( .01 )*FI l<-( Al+Bl*V-x-V) 

G4=(-G1*(2.^-D1*(W-W0) )+G2*(2.*D1*(W-W0) ) )/(G2*G2) 

G5= (C1+D1*( W-WO)* (W-WO ) )/ ( l.+Dl*(W-WO)*(W-WO) ) 

G7=l.+W/W0 

Pl = ( ( l.+Dl^<-( W-WO) -X (W-WO) )x- 2,*D1+4.*D1*( W-WO) *D1*{ W-WO) )/(G2*G2) 
P2=(-(8.*Dl-x-(W-WO)*(W-WO) )*D1) /(G2*G2) 

P3= (-( C 1 +D 1-x- ( W-WO )^-( W-WO) )*2,*D1-4.*D1*(W-W0)*D1*( W-WO) )/ (G2*G2) 
P4= ( G 1*8. *D1* ( W-WO) *D1*( W-WO) ) /(G2*G2*G2) 

P9=FI2*( ( (-G7*( .1/W0)+G6*( l./WO) )*2.*(1./W0) )/(G7*G7*G7) ) 
HV=XLAM*C0SU+XMU*S!NU+( .01 )*FI1*( 2.*B1*V)*G5*( l.-.25*COSUMQ) 

HU = -XLA,M*V*SINU+XMU*V*COSU+G3*G5*( .25*SINUMO ) 

HVV=( .01)*2.*B1*FI 1*G5*( 1.-.25*C0SUMQ) 

HVU=-XLAM*SINU+XMU*COSU+G5*.Ol*FIl*2.*Bl*V*( .25*SINUMO) 
hUU=-XLAM*V*COSU-XMU*V*SINU+G3*G5*( .25*COSUMO) 

HW=G3*( l.-.25*( COSUMQ) )*G4+(Fl2*( (G7*(.l/WO) )-G6*(l./WO) ) )/(G7*G7) 

HWW= G3*( ( l.-.25*COSUMQ) *’(P1 + P2+P3+P4) )+P9 

D£N=HVU*HVU-HUU*HVV 

GUM1=-HV*HVU+HU*HVV 

GUM2=HVU*( -HU ) +HV*HUU 

DO 981 1=1,7 

DU=GUM1/DEN 

DV=GUM2/DEN 

U=U-i-DU 

V=V+DV 

COSU=COSF(U ) 

SINU=SINF(U) 
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COSUMQ=COSF (U-Q) 

SINUMQ=SINF (U-Q ) 

G1=C1+D1*(W-W0) *( W-WO) 

G2 = l.+Dl*(W-WO)-«-(W-WO) 

G3=( .01)-«-FI 1«-(A1+B1*V^^V) 

G4=(-G1*(2.*D1^-(W-W0) )+G2*(2.*D1*( W-WO) ) )/(G2*G2) 
G5=(Cl+Dl^HW-WO)*(W-WO) )/( 1 . +D 1* ( W-WO ) * ( W-WO ) ) 
F1V=(A1+B1*V^:-V)*( .01 ) «• { 1 . - . 2 5*COSUMQ ) *G 1/G2 
F2V=(l.+.l*W/WO)/( l.+W/WO) 

HV = XLAM*COSU+XMU*SINU+( .01 ) *F I 1* C 2 . *B1*V ) *G5* ( 1 . - . 2 5 *COSUMO ) 
HU = -XLAM«-V*S IMU+XMU*V*COSU+G3*G5* ( .2 5*5 I NUMO ) 

HVV=( .01 )*2.*-Bl*FI 1*G5*( 1 . - . 2 5*COSUMQ ) 

HVU = -XLAM*SINU+XMU*COSU+G5*.Oi*FI1*2.*B1*V*( .2 5*5 1 NUMO) 
HUU=-XLAM*V*C05U-XMU*V*5INU+G3*G5*( .25*C05UMQ) 
DEN=HVU*HVU-HUU*HVV 
GUM1=-HV*HVU+HU*HVV 
GUM2=HVU* ( -HU ) +HV*HUU 
981 CONTINUE 
RETURN 
END 
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Appendix II 



In this appendix additional information about the 
paths described in section 7 and the path discussed in 
section 8 x^ill be given. For the first three paths 
the position (xjy) and the controls u,v,w for the sub- 
raarine v:ill be given at twelve -hour intervals. This 
information I'n.ll be given for Path IV and in addition 
the submarine’s location and controls will be included 
for the time steps immediately preceeding and following 
the corner. 



Path I 



Time 


X 


Z 


OoO 


0.0 


0.0 


12.0 


104.5 


- 7^.6 


24.0 


214.8 


- 140,4 


36.0 


330.0 


- 197.2 


48.0 


449.2 


- 244.8 


60.0 


571.5 


- 283.1 


72.0 


696,1 


- 312.5 


84.0 


822,3 


- 333oO 


96.0 


9^9.3 


- 3^5.1 



u 


V 


w 


0.7 


10.7 


207.8 


0.6 


10.7 


207.8 


0.5 


10.7 


207.8 


0.4 


10.7 


207.8 


0.3 


10.7 


207.7 


0.3 


10.7 


207.5 


0.2 


10.7 


207.3 


0.1 


10.6 


207.1 


0.1 


10.6 


206.8 
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Time 


X 


Z 


u 


V 


w 


108.0 


1076.5 


- 3^9.2 


0.0 


10.6 


206.5 


120.0 


1203.5 


- 3^5.7 


0.1 


10.6 


206.2 


132.0 


1329.8 


- 335.1 


0.1 


10.5 


205.9 


144.0 


1455.0 


- 318.0 


0.2 


10.5 


205.5 


156.0 


1578.9 


- 294.8 


0.2 


10.5 


205.2 


168.0 


1701.3 


- 266.0 


0.3 


10.5 


204.9 


180.0 


1822.0 


- 232.1 


0.3 


10.4 


204.6 


192.0 


1941.0 


- 193.6 


0.3 . 


10.4 , 


, 204.3 


204.0 


2058.2 


- 150.8 


0.4 


10.4 


204.0 


216.0 


2173.6 


- 104.2 


0.4 


10.4 


203.7 


228.0 


2287.2 


- 5^.1 


0.4 


10.3 


203.5 


240.0 


2400.0 


0.0 


0.5 


11.1 


203.7 






Path II 






Time 


X 


Z 


u 


V 


w 


0.0 


0.0 


0.0 


- 0.6 


10.1 


1000.0 


12.0 


102.4 


- 66.1 


- 0.5 


10.2 


1000.0 


24.0 


210.1 


- 124.3 


- 0.5 


10.2 


1000.0 


36.0 


322.4 


- 174.3 


- 0.4 


10.3 


1000.0 


48.0 


438.6 


- 216.3 


•- 0.3 


10.3 


1000,0 


60.0 


558.0 


- 250.4 


- 0.2 


10.4 


1000.0 


72.0 


680.0 


- 276.7 


- 0.2 


10.4 


1000.0 
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T3Jtl3 


X 


Z 


u 


V 


w 


aij'.o 


804.0 


- 295.6 


- 0.1 


10.5 


1000.0 


96.0 


929.4 


- 307.^ 


- 0.1 


10.5 


1000.0 


1C6.0 


1055.7 


- 312.3 


0.0 


10.5 


1000.0 


120.0 


1182.4 


- 310.9 


0.0 


10.6 


1000.0 


132.0 


1309.2 


- 303.7 


0.1 


10.6 


1000.0 


li|4.0 


1^'35.7 


- 290.9 


0.1 


10.6 


1000.0 


156.0 


1561.7 


- 273.2 


0.2 


10.6 


1000.0 


168.0 


1686.9 


- 250.8 


0.2 


10.6 


1000.0 


180.0 


1811.2 


- 224.1 


0.2 


10.6 


1000.0 


192.0 


193^.5 


- 193.7 


0.3 


10.6 


1000.0 


204.0 


2056.6 


- 159.8 


0.3 


10.6 


1000.0 


216.0 


2177.6 


- 122.7 


0.3 


10.5 


1000.0 


228.0 


2290.2 


- 64.0 


0.5 


10.6 


705.2 


240.0 


2400.0 


0.0 


0.5 


10.7 


650.1 






Path 


III 






Time 


X 


Z 


u 


V 


w 


0.0 


0.0 


0.0 


-'0.7 


10.8 


527.4 


12.0 


106.1 


- 75.^ 


- 0.6 


10.9 


526.4 


24.0 


218.1 


- 141.9 


- 0.5 


10.9 


526.6 


36.0 


33^.9 


- 199.1 


- 0.4 


10.8 


528.1 


48.0 


455.8 


- 2k6,9 


- 0.3 


10.8 


530.7 
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Tirao 


X 


1 


u 


V 


W 


60.0 


579.6 


- 285.4 


- 0.3 


10.8 


53^.6 


72.0 


705.5 


- 31^.6 


- 0.2 


10.8 


539.6 


84.0 


832.8 


- 33^.8 


- 0.1 


10.7 


545.9 


96.0 


960.6 


- 346.5 


- 0.1 


10.7 


553.3 


108.0 


1088.3 


- 350.1 


0.0 


10.6 


561.9 


120.0 


1215.4 


- 3^^6.1 


0.1 


10.6 


571.7 


132.0 


1341.5 


- 335.0 


0.1 


■ 10.5 


582.7 


144,0 


1466.3 


- 317.4 


0.2 


10.5 


594.9 


156.0 


1589.5 


- 293.9 


0.2 


10.4 


608.3 


168.0 


1710.9 


- 264.8 


0.3 


10.4 


622.9 


180.0 


1830.5 


- 230.7 


0.3 


10.3 


638.7 


192.0 


1948.1 


- 192.1 


0.3 


10.3 


655.8 


204.0 


2063.8 


- 149.4 


0.4 


10.3 


674.2 


216.0 


2177.6 


- 103.0 


0.4 


10.2 


693.9 


228.0 


2289.6 


- 53.2 


0.4 


10.2 


715.0 


240.0 


2400,1 


0.0 


0.5 


10.4 


677.9 






Path IV 






Ties 


X 


y 


u 


V ■ 


w 


0.0 


0.0 


0.0 


•- 1.1 


11.0 


200.8 


• 12.0 


70.2 


- 111.9 


- 1.0 


11.0 


200.8 
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Time 



as Z 

24.0 149.6 - 218.0 

36.0 237.7 - 317.5 

48.0 334.1 - 409.4 

60.0 437.9 - 493.1 

72.0 548.4 - 568.0 

84.0 664.6 - 633.5 

96.0 785.7 - 689.6 

108.0 910.5 - 736.0 

120.0 1038.3 - 772.7 

132.0 1168.3 - 799.8 

144.0 1302.0 - 816.3 

156.0 1453.5 - 811.1 

158.0 1482.9 - 804.7 

160.0. 1511.5 - 796.0 

162.0 1539.4 - 784.8 

164.0 1565.1 - 771.6 

166.0 1591.5 - 757.2 

168.0 1617.3 - 742.0 

180.0 1759.6 - 632.3 

192.0 1889.7 - 507.8 

204.0 2017.5 - 381.1 

216.0 2145.1 - 254.1 

228.0 2272.6 - 127.0 

240.0 2400.0 0.1 



u 


V 


w 


0.9 


11.1 


200.8 


0.8 


11.1 


200.8 


0.7 


11.1 


200.8 


0.6 


11.1 


200.8 


0.6 


11.1 


200.9 


0.5 


11.1 


200.8 


0.4 


11.1 


200.8 


0.3 


11.1 


200.8 


0.2 


11.1 


200.8 


0.2 


11.1 


200.8 


0.1 


11.5 


201.1 


0.2: 


- 15.0 


202.9 


0.3 


15.0 


204.5 


0.3 


15.0 


207.0 


0.4 


15.0 


211.4 


0.5 


15.0 


1000.0 


0.5 


15.0 


1000.0 


0.6 


15.0 


1000.0 


0.7 


15.0 


1000.0 


0.8 


15.0 


1000.0 


0.8 


15.0 


1000.0 


0.8 


15.0 


1000.0 


0.8 


15.0 


1000.0 


0.8 


15.0 


1000.0 
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